module Read_Parameters
!
! Module for reading parameter files
!
! Read the file with:
!
! CALL read_parameterfile(filename [,err_unit=...] [,verbose=...]
!                         [,output_parameterfile=...],
!                         [, restarting =...] [, restartfile=...] )
!
! And extract parameters with:
!
! CALL get_parameter(parametername,variable [,default=...] [,nelements=...]
!                    [,required=...] [,allowed =...] [, range=...]
!                    [, is_set=...],  )
!
! which reads a parameter into a variable. If the parametername is not found
! in the parameter file the default value is used. If no default is
! specified and the parameter is not in the file an error occurs.
!
! Lines that cannot be translated into the same type as the
! variable will cause the program to terminate with an error message.
!
! 'filename' and 'parametername' are strings. 'variable' can be a scalar or
! one dimensional array of type real, integer, logical, double precision
! or string. err_unit is the unit number for error/diagnostic messages
! and setting verbose=.TRUE. outputs parameter values to this unit. If
! err_unit is not specified stdout is used.
!
! Do not use enclosing quotes for string parameters unless the quotes are
! part of the string.
!
! The optional parameter 'nelements' returns the number of elements in an array
! parameter (see below).
!
! If required=.FALSE. the parameter is not read - this is intended for the
! case where the setting of one parameter may make another unnecessary.
! With required=false no error occurs if the parameter is not found
! or if the parameter is present but not required.
!
! If required=true but a default is specified, then no error will occur if
! the parameter is not present in the parameterfile. The default setting
! is required=.true. This is really intended for cases where one parameter
! setting means that another parameter is unnecessary.
!
! 'range' (a 2 element array) lets you specify maximum and minimum acceptable
! values for the parameter and can be used for real/double/integer parameters.
! e.g. use range=(/ 0, 10 /) to specify that values 0-10 inclusive are ok.
!
! 'allowed' (an array) lets you specify which values are acceptable for
! integer parameters, e.g. allowed= (/ 0, 5, 2, 9 /)
!
! 'is_set' (a logical) returns true if the parameter is explicitly set in
! the parameter file and false if its not (eg if a default value is used
! or the parameter is missing and has required=.false.)
!
! Once all the parameters have been read, do
!
! CALL finalize_parameters( [abort_on_unused=...] [,abort_on_default] )
!
! to check for unused parameters in the file and generate warnings about
! any default values that have been used. If abort_on_unused is true then
! unused parameters in the file will cause an error. This is to try to catch
! typos in the names of optional parameters.
!
! If abort_on_default is true then the program will exit if any parameters
! (other than those with required=false) are not explicitly specified in the
! parameter file.
!
! The default is abort_on_default = abort_on_unused = .FALSE.
!
! Parameter files must consist of lines with
!
! <parameter_name> = <parameter_value>
!
! Blank lines are ignored and comments can be included by preceding them with
! a hash.
!
! Array parameters can be included. For example, for a 5 element
! array of integers called 'iarr' the line in the parameter file would be:
!
! iarr = [1, 5, 8, 12, 9]
!
! To specify a default value for this array parameter you could include
! something like "default=(/ 5,6,7,8,9 /)" in the call to get_parameter.
!
! Arrays must be in square brackets with commas between elements.
! Note that its not possible to read an array of strings if some of
! the strings contain commas, because the comma is used as a separator.
!
! A backslash (\) as the last character on a line indicates a
! continuation line.
!
! Parameter names are case insensitive.
!
  use log
  use get_unitmod
  implicit none
  save
  private

  public :: read_parameterfile
  public :: get_parameter
  public :: finalize_parameters

  ! Array to store lines from file
  ! ncharmax needs to be quite large if we have string array parameters
  ! and/or continuation lines.
  integer :: nlines = -1
  integer, parameter :: nlinemax = 256, ncharmax = 4000
  character(LEN=ncharmax), dimension(nlinemax) :: line

  ! Arrays used to generate warnings/errors
  integer :: ndefault
  character(LEN=ncharmax), dimension(nlinemax) :: default_param
  logical, dimension(nlinemax) :: used_param

  ! Unit number for errors 
  integer :: err_params = -1

  character(LEN=2), parameter :: cont_line = "\ "
  
  ! Generic routine to extract parameters
  interface get_parameter
     ! These shouldn't be called directly - use 'CALL get_parameter(...)'
     module procedure get_integer_parameter
     module procedure get_logical_parameter
     module procedure get_real_parameter
     module procedure get_string_parameter
     module procedure get_double_parameter
     module procedure get_integer_array_parameter
     module procedure get_real_array_parameter
     module procedure get_double_array_parameter
     module procedure get_logical_array_parameter
     module procedure get_string_array_parameter
  end interface

  logical :: verbose_params = .true.

  ! New parameterfile to write
  integer :: out_lun, restart_lun
  character(LEN=500) :: output_filename
  logical :: write_file = .false.

  ! Restart files
  logical :: make_restartfile = .false. , read_restartfile = .false.

  contains

    subroutine read_parameterfile(fname,err_unit,verbose,output_parameterfile, &
         restartfile, restarting)
!
! This reads in the parameter file and checks for invalid lines.
! It also throws away any comments and blank lines and stitches together
! continuation lines.
!
      use log
      implicit none
      character(LEN=*), intent(IN) :: fname
      character(LEN=20) :: fstr
      integer :: ios, i
      character(LEN=ncharmax) :: str, nextline
      integer, optional :: err_unit
      logical, optional :: verbose, restarting
      logical :: failed
      integer :: lun_params
      character(LEN=*), optional :: output_parameterfile, restartfile
      logical :: cont

      used_param = .false.
      ndefault = 0

      if(present(verbose))verbose_params=verbose
      if(present(err_unit))err_params=err_unit

      ! Check if restarting
      if(present(restartfile))then
         restart_lun = get_unit()
         if(.not.present(restarting)) call readparams_abort( &
              'read_parameterfile() : Keyword "restartfile" was specified without keyword "restarting"')
         if(restarting)then
            make_restartfile = .false.
            read_restartfile = .true.
            open(unit=restart_lun, file=restartfile, status='old', form='unformatted',action='read')
         else
            make_restartfile = .true.
            read_restartfile = .false.
            open(unit=restart_lun, file=restartfile, status='unknown',form='unformatted',action='write')
         end if
      else
         make_restartfile = .false.
         read_restartfile = .false.
      end if
      if(present(restarting).and.(.not.present(restartfile))) call readparams_abort( &
           'read_parameterfile() : Keyword "restarting" was specified without keyword "restartfile"')


      ! Open the output file if necessary
      if(present(output_parameterfile))then
         out_lun = get_unit()
         open(unit=out_lun,file=output_parameterfile,status='unknown', &
              form='formatted')
         write_file = .true.
      else
         write_file=.false.
      end if

      ! Make format string to read a line
      write(fstr,'("(a",i4,")")')ncharmax

      ! Find a free unit number
      lun_params = get_unit()

      ! Read lines until end of file
      nlines = 0
      ios = 0
      open(unit=lun_params,form='formatted',status='old',file=fname,iostat=ios)
      if(ios.ne.0)call readparams_abort( &
           "Unable to open parameter file: "//trim(fname))
      do while(ios.eq.0)
         read(unit=lun_params,fmt=fstr,iostat=ios)line(nlines+1)
         if(ios.eq.0)then
            nlines=nlines+1
            ! Check if it starts with a # (comment line)
            str=trim(adjustl(line(nlines)))
            if(index(str,"#").eq.1)then
               nlines=nlines-1
               cycle
            else
               if(len_trim(line(nlines)).eq.0)then
                  nlines=nlines-1
                  cycle
               end if
            endif
            ! Check for continuation lines
            if(str(len_trim(str):len_trim(str)).eq.cont_line(1:1))then
               cont = .true.
               do while(cont)
                  read(unit=lun_params,fmt=fstr,iostat=ios)nextline
                  if(ios.ne.0)call readparams_abort( &
                       'Unable to read continuation line')
                  if(index(nextline,cont_line(1:1)).gt.0)then
                     cont=.true.
                  else
                     cont=.false.
                  end if
                  str = line(nlines)

                  if(len_trim(str(1:len_trim(str)-1)//trim(adjustl(nextline))).gt.ncharmax) &
                       call readparams_abort('Line too long after continuation - increase ncharmax!')

                  line(nlines) = &
                       str(1:len_trim(str)-1)//trim(adjustl(nextline))
               end do
            end if
         end if
         if(nlines.eq.nlinemax) &
              call readparams_abort('Too many lines in parameter file!')
      end do
      close(lun_params)

      if(nlines.lt.1)&
           call readparams_abort('No lines read from parameter file!')

      failed = .false.
      do i=1,nlines,1
         if(index(line(i),"=").eq.0)then
            if(err_params.eq.-1)then
               write(logunit,*)'Can''t interpret line: ',trim(adjustl(line(i)))
            else
               write(err_params,*) &
                    'Can''t interpret line: ',trim(adjustl(line(i)))
            endif
            call readparams_abort('Failed reading file: '//trim(fname))
         end if
      end do

      if(verbose_params)then
         if(err_params.eq.-1)then
            write(logunit,*)'  -read ',nlines,' lines from the parameter file ' &
                 //trim(fname)
         else
            write(err_params,*) &
                 '  -read ',nlines,' lines from the parameter file '&
                 //trim(fname)
         end if
      end if

      return

    end subroutine read_parameterfile

!
! The following routines get a parameter from the lines read into memory
! or assign a default value to the parameter.
!
! They should not be called directly - use the generic interface:
!
! CALL get_parameter(name,param,default=...)
!
! which will automatically select the appropriate routine for the parameter
! specified.
!
    subroutine get_integer_parameter(name,iparam,default,required, &
                                     range, allowed, is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      integer, intent(OUT) :: iparam
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios, itemp
      integer, optional :: default
      logical, optional :: required
      integer, dimension(2), optional :: range
      integer, dimension(:), optional :: allowed
      logical, optional :: is_set

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 ' unable to find parameter '//trim(name)//' in parameter file')
         else
            iparam = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.ge.0)then
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         parval = str(ifirst:ilast)

         read(parval,*,iostat=ios)iparam
         if(present(is_set)) is_set = .true.

         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))
         
         if(verbose_params.and.err_params.eq.-1)write(logunit,*)'  -parameter ',trim(name),' = ',iparam
         if(verbose_params.and.err_params.ne.-1)write(err_params,*)'  -parameter ',trim(name),' = ',iparam

         used_param(jline) = .true.
      end if

      if(present(range))then
         if(iparam.lt.range(1).or.iparam.gt.range(2)) call readparams_abort( &
              '  -parameter outside allowed range: '//trim(name))
      end if

      if(present(allowed))then
         if(.not.any(allowed.eq.iparam)) call readparams_abort( &
              '  -parameter has illegal value: '//trim(name))
      end if

      if(write_file)then
         write(out_lun,'(a30,a3,i13)')trim(name)," = ",iparam
      end if

      if(make_restartfile) write(restart_lun)iparam
      if(read_restartfile) then
         read(restart_lun)itemp
         if(itemp.ne.iparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if
      
      return
    end subroutine get_integer_parameter


    subroutine get_logical_parameter(name,lparam,default,required,is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      logical, intent(OUT) :: lparam
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      logical, optional :: default
      logical, optional :: required
      logical, optional :: is_set
      logical :: ltemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort(' parameter file has not been read!')

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(default))then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            lparam = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      endif

      if(jline.gt.0)then
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         parval = str(ifirst:ilast)

         read(parval,*,iostat=ios)lparam
         if(present(is_set)) is_set = .true.

         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))
         
         if(verbose_params.and.err_params.eq.-1)write(logunit,*)'  -parameter ',trim(name),' = ',lparam
         if(verbose_params.and.err_params.ne.-1)write(err_params,*)'  -parameter ',trim(name),' = ',lparam
         
         used_param(jline) = .true.
      end if

      if(write_file)then
         write(out_lun,'(a30,a3,l13)')trim(name)," = ",lparam
      end if

      if(make_restartfile) write(restart_lun)lparam
      if(read_restartfile) then
         read(restart_lun)ltemp
         if(ltemp.neqv.lparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if

      return
    end subroutine get_logical_parameter


    subroutine get_real_parameter(name,rparam,default,required,range,is_set)

      use log
      implicit none
      real :: rtemp
      character(LEN=*), intent(IN) :: name
      real, intent(OUT) :: rparam
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      real, optional :: default
      logical, optional :: required, is_set
      real, dimension(2), optional :: range
      
      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(default))then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            rparam = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         endif
      end if

      if(jline.gt.0)then
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         parval = str(ifirst:ilast)

         read(parval,*,iostat=ios)rparam
         if(present(is_set)) is_set = .true.

         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))

         if(verbose_params.and.err_params.eq.-1)write(logunit,*)'  -parameter ',trim(name),' = ',rparam
         if(verbose_params.and.err_params.ne.-1)write(err_params,*)'  -parameter ',trim(name),' = ',rparam
         
         used_param(jline) = .true.
      end if

      if(present(range))then
         if(rparam.lt.range(1).or.rparam.gt.range(2)) call readparams_abort( &
              '  -parameter outside allowed range: '//trim(name))
      end if

      if(write_file)then
         write(out_lun,'(a30,a3,es13.5)')trim(name)," = ",rparam
      end if

      if(make_restartfile) write(restart_lun)rparam
      if(read_restartfile) then
         read(restart_lun)rtemp
         if(rtemp.ne.rparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if

      return
    end subroutine get_real_parameter


    subroutine get_double_parameter(name,dparam,default,required,range,is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      double precision, intent(OUT) :: dparam
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      double precision, optional :: default
      logical, optional :: required, is_set
      double precision, dimension(:), optional :: range
      double precision :: dtemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(default))then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            dparam = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         endif
      end if

      if(jline.gt.0)then
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         parval = str(ifirst:ilast)

         read(parval,*,iostat=ios)dparam
         if(present(is_set)) is_set = .true.

         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))

         if(verbose_params.and.err_params.eq.-1)write(logunit,*)'  -parameter ',trim(name),' = ',dparam
         if(verbose_params.and.err_params.ne.-1)write(err_params,*)'  -parameter ',trim(name),' = ',dparam
         
         used_param(jline) = .true.
      end if

      if(present(range))then
         if(dparam.lt.range(1).or.dparam.gt.range(2)) call readparams_abort( &
              '  -parameter outside allowed range: '//trim(name))
      end if

      if(write_file)then
         write(out_lun,'(a30,a3,d13.5)')trim(name)," = ",dparam
      end if

      if(make_restartfile) write(restart_lun)dparam
      if(read_restartfile) then
         read(restart_lun)dtemp
         if(dtemp.ne.dparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if

      return
    end subroutine get_double_parameter


    subroutine get_string_parameter(name,sparam,default,required,is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      character(LEN=*), intent(OUT) :: sparam
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast
      character(LEN=*), optional :: default
      logical, optional :: required, is_set
      character(LEN=ncharmax) :: stemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(default))then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            sparam=default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",trim(default)
         end if
      end if

      if(jline.gt.0)then
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         parval = str(ifirst:ilast)

         sparam = trim(adjustl(parval))
         if(present(is_set)) is_set = .true.

         if(verbose_params.and.err_params.eq.-1)write(logunit,*)'  -parameter ',trim(name),' = ',trim(sparam)
         if(verbose_params.and.err_params.ne.-1)write(err_params,*)'  -parameter ',trim(name),' = ',trim(sparam)
         
         used_param(jline) = .true.
      end if
         
      if(write_file)then
         write(out_lun,'(a30,a3,a)')trim(name)," = ",trim(sparam)
      end if

      stemp = trim(sparam)
      if(make_restartfile) write(restart_lun)stemp
      if(read_restartfile) then
         read(restart_lun)stemp
         if(trim(stemp).ne.trim(sparam))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if

      return
    end subroutine get_string_parameter

!
! This stops the program in the event of an error
!
    subroutine readparams_abort(str)

      use log
      implicit none
      character(LEN=*) :: str

      if(err_params.eq.-1)then
         write(logunit,*)'  -read_parameters(): FATAL'
         write(logunit,*)'  '//trim(str)
         call flush(0)
         call flush(5)
         call flush(6)
         stop
      else
         write(err_params,*)'  -read_paramaters(): FATAL'
         write(err_params,*)'  '//trim(str)
         call flush(err_params)
         stop
      endif

    end subroutine readparams_abort

!
! This checks that all the lines have been used and warns about default
! values if verbose=true.
!
    subroutine finalize_parameters(abort_on_unused, abort_on_default)

      use log
      implicit none
      logical, optional :: abort_on_unused, abort_on_default
      integer :: i
      integer :: nunused

      if(verbose_params)then
         if(err_params.ne.-1)then
            write(err_params,*)'  -checking parameters...'
         else
            write(logunit,*)'  -checking parameters...'
         end if

         do i=1,ndefault,1
            if(err_params.ne.-1)then
               write(err_params,*)"* Parameter set to default: ", &
                    trim(default_param(i))
            else
               write(logunit,*)"* Parameter set to default: ", &
                    trim(default_param(i))
            end if
         end do
         
         if(ndefault.gt.0)then
            if(err_params.ne.-1)then
               !write(err_params,*)" "
            else
               !write(logunit,*)" "
            end if
         end if
      end if

      nunused = 0
      do i=1,nlines,1
         if(.not.used_param(i))then
            nunused = nunused + 1
            if(verbose_params)then
               if(err_params.ne.-1)then
                  write(err_params,*)"WARNING - Unused parameter: ", &
                       trim(line(i))
               else
                  write(logunit,*)"WARNING - Unused parameter: ", &
                       trim(line(i))
               end if
            end if
         end if
      end do

      if(nunused.gt.0.and.verbose_params)then
         if(err_params.ne.-1)then
            write(err_params,*)" "
         else
            write(logunit,*)" "
         end if
      end if

      if(verbose_params)then
         if(err_params.ne.-1)then
            write(err_params,*)"  -there are ",nunused, &
                 " unused parameter(s) in the file"
            write(err_params,*)"  -there are ",ndefault, &
                 " parameter(s) set to default values"
         else
            write(logunit,*)"  -there are ",nunused, &
                 " unused parameter(s) in the file"
            write(logunit,*)"  -there are ",ndefault, &
                 " parameter(s) set to default values"
         end if
      end if

      if(present(abort_on_default))then
         if(abort_on_default.and.ndefault.gt.0)then
            call readparams_abort( &
                 "  -there are parameters with no entry in the parameter "// &
                 "file and abort_on_default=.true.")
         end if
      end if

      if(present(abort_on_unused))then
         if(abort_on_unused.and.nunused.gt.0)then
            call readparams_abort("  -there are unused parameters!")
         end if
      end if

      if(verbose_params)then
         if(err_params.ne.-1)then
            write(err_params,*)"  -finished checking parameters"      
         else
            write(logunit,*)"  -finished checking parameters"      
         end if
      end if

      ! Close files
      if(write_file)close(out_lun)
      if(make_restartfile.or.read_restartfile)close(restart_lun)

      return
    end subroutine finalize_parameters

!
! The following routines are to deal with the case where the parameter is
! a 1D array
!
    subroutine get_integer_array_parameter(name,array,default,nelements, &
         required, is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      integer, intent(OUT), dimension(:) :: array
      integer, optional, intent(OUT) :: nelements
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      integer, optional, dimension(:) :: default
      integer :: ilower, iupper, narray, nparam
      integer :: i
      character(LEN=256) :: fmtstring
      logical, optional :: required, is_set
      integer :: ntemp
      integer, allocatable, dimension(:) :: atemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      ilower = lbound(array,1)
      iupper = ubound(array,1)
      narray = iupper-ilower+1

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            array = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.gt.0)then

         if(present(is_set)) is_set = .true.

         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         ! Get the parameter value minus any extra spaces at either end
         parval = trim(adjustl(str(ifirst:ilast)))
      
         ! Check that it is enclosed in square brackets
         if(parval(1:1).ne."[".or. &
              parval(len_trim(parval):len_trim(parval)).ne."]") &
              call readparams_abort('Arrays must be in square brackets'// &
              ' with elements separated by commas')

         str = trim(adjustl(parval(2:len_trim(parval)-1)))

         ! Count how many elements there are
         nparam = 1
         do i=1,len_trim(str),1
            if(str(i:i).eq.",")nparam=nparam+1
         end do

         if(nparam.gt.narray) &
              call readparams_abort("Array parameter "//trim(name) &
              //" does not fit in the supplied array!")
         
         read(str,*,iostat=ios)array(ilower:ilower+nparam-1)
         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))

         if(verbose_params.and.err_params.eq.-1) &
              write(logunit,*)'  -parameter ',trim(name),' = ',array
         if(verbose_params.and.err_params.ne.-1) &
              write(err_params,*)'  -parameter ',trim(name),' = ',array

         used_param(jline) = .true.

      else
         nparam = ubound(default,1)-lbound(default,1)+1
      end if

      if(present(nelements))nelements = nparam

      if(write_file)then
         if(nparam.gt.1)then
            write(fmtstring, &
                 '("(a30,"" = "",a1,""["",",i13,"(i13,"",""),i13,""]"")")') &
                 nparam-1
         else
            fmtstring = '(a30,a3,"[ ",i13," ]")'
         end if
         write(out_lun,fmtstring)trim(name)," = ",array(1:nparam)
      end if

      if(make_restartfile)then
         write(restart_lun)nparam
         write(restart_lun)array(1:nparam)
      end if
      if(read_restartfile) then
         read(restart_lun)ntemp
         if(ntemp.ne.nparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         allocate(atemp(1:nparam))
         read(restart_lun)atemp(1:nparam)
         if(any(atemp(1:ntemp).ne.array(1:nparam)))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         deallocate(atemp)
      end if

      return
    end subroutine get_integer_array_parameter


    subroutine get_real_array_parameter(name,array,default,nelements, &
         required,is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      real, intent(OUT), dimension(:) :: array
      integer, optional, intent(OUT) :: nelements
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      real, optional, dimension(:) :: default
      integer :: ilower, iupper, narray, nparam
      integer :: i
      character(LEN=256) :: fmtstring
      logical, optional :: required, is_set
      integer :: ntemp
      real, allocatable, dimension(:) :: atemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      ilower = lbound(array,1)
      iupper = ubound(array,1)
      narray = iupper-ilower+1

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            array = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.gt.0)then

         if(present(is_set)) is_set = .true.

         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         ! Get the parameter value minus any extra spaces at either end
         parval = trim(adjustl(str(ifirst:ilast)))
      
         ! Check that it is enclosed in square brackets
         if(parval(1:1).ne."[".or. &
              parval(len_trim(parval):len_trim(parval)).ne."]") &
              call readparams_abort('Arrays must be in square brackets'// &
              ' with elements separated by commas')
         
         str = trim(adjustl(parval(2:len_trim(parval)-1)))

         ! Count how many elements there are
         nparam = 1
         do i=1,len_trim(str),1
            if(str(i:i).eq.",")nparam=nparam+1
         end do
         
         if(nparam.gt.narray) &
              call readparams_abort("Array parameter "//trim(name) &
              //" does not fit in the supplied array!")
         
         read(str,*,iostat=ios)array(ilower:ilower+nparam-1)
         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))
         
         if(verbose_params.and.err_params.eq.-1) &
              write(logunit,*)'  -parameter ',trim(name),' = ',array
         if(verbose_params.and.err_params.ne.-1) &
              write(err_params,*)'  -parameter ',trim(name),' = ',array

         used_param(jline) = .true.
      else
         nparam = ubound(default,1)-lbound(default,1)+1
      end if

      if(present(nelements))nelements = nparam

      if(write_file)then
         if(nparam.gt.1)then
            write(fmtstring, &
                 '("(a30,"" = "",a1,""["",",i13,"(es13.5,"",""),es13.5,""]"")")') &
                 nparam-1
         else
            fmtstring = '(a30,a3,"[ ",es13.5," ]")'
         end if
         write(out_lun,fmtstring)trim(name)," = ",array(1:nparam)
      end if

      if(make_restartfile)then
         write(restart_lun)nparam
         write(restart_lun)array(1:nparam)
      end if
      if(read_restartfile) then
         read(restart_lun)ntemp
         if(ntemp.ne.nparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         allocate(atemp(1:nparam))
         read(restart_lun)atemp(1:nparam)
         if(any(atemp(1:ntemp).ne.array(1:nparam)))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         deallocate(atemp)
      end if

      return
    end subroutine get_real_array_parameter

    subroutine get_double_array_parameter(name,array,default,nelements, &
         required,is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      double precision, intent(OUT), dimension(:) :: array
      integer, optional, intent(OUT) :: nelements
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      double precision, optional, dimension(:) :: default
      integer :: ilower, iupper, narray, nparam
      integer :: i
      character(LEN=256) :: fmtstring
      logical, optional :: required, is_set
      integer :: ntemp
      double precision, allocatable, dimension(:) :: atemp

      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      ilower = lbound(array,1)
      iupper = ubound(array,1)
      narray = iupper-ilower+1

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            array = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.gt.0)then

         if(present(is_set)) is_set = .true.

         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)
         
         ! Get the parameter value minus any extra spaces at either end
         parval = trim(adjustl(str(ifirst:ilast)))
      
         ! Check that it is enclosed in square brackets
         if(parval(1:1).ne."[".or. &
              parval(len_trim(parval):len_trim(parval)).ne."]") &
              call readparams_abort('Arrays must be in square brackets'// &
              ' with elements separated by commas')

         str = trim(adjustl(parval(2:len_trim(parval)-1)))

         ! Count how many elements there are
         nparam = 1
         do i=1,len_trim(str),1
            if(str(i:i).eq.",")nparam=nparam+1
         end do

         if(nparam.gt.narray) &
              call readparams_abort("Array parameter "//trim(name) &
              //" does not fit in the supplied array!")

         read(str,*,iostat=ios)array(ilower:ilower+nparam-1)
         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))

         if(verbose_params.and.err_params.eq.-1) &
              write(logunit,*)'  -parameter ',trim(name),' = ',array
         if(verbose_params.and.err_params.ne.-1) &
              write(err_params,*)'  -parameter ',trim(name),' = ',array

         used_param(jline) = .true.
      else
         nparam = ubound(default,1)-lbound(default,1)+1
      end if

      if(present(nelements))nelements = nparam

      if(write_file)then
         if(nparam.gt.1)then
            write(fmtstring, &
                 '("(a30,"" = "",a1,""["",",i13,"(d13.5,"",""),d13.5,""]"")")') &
                 nparam-1
         else
            fmtstring = '(a30,a3,"[ ",d13.5," ]")'
         end if
         write(out_lun,fmtstring)trim(name)," = ",array(1:nparam)
      end if

      if(make_restartfile)then
         write(restart_lun)nparam
         write(restart_lun)array(1:nparam)
      end if
      if(read_restartfile) then
         read(restart_lun)ntemp
         if(ntemp.ne.nparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         allocate(atemp(1:nparam))
         read(restart_lun)atemp(1:nparam)
         if(any(atemp(1:ntemp).ne.array(1:nparam)))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         deallocate(atemp)
      end if

      return
    end subroutine get_double_array_parameter


    subroutine get_logical_array_parameter(name,array,default,nelements, &
         required, is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      logical, intent(OUT), dimension(:) :: array
      integer, optional, intent(OUT) :: nelements
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      logical, optional, dimension(:) :: default
      integer :: ilower, iupper, narray, nparam
      integer :: i
      character(LEN=256) :: fmtstring
      logical, optional :: required, is_set
      integer :: ntemp
      logical, allocatable, dimension(:) :: atemp
   
      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      ilower = lbound(array,1)
      iupper = ubound(array,1)
      narray = iupper-ilower+1

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            array = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.gt.0)then

         if(present(is_set)) is_set = .true.
         
         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         ! Get the parameter value minus any extra spaces at either end
         parval = trim(adjustl(str(ifirst:ilast)))
      
         ! Check that it is enclosed in square brackets
         if(parval(1:1).ne."[".or. &
              parval(len_trim(parval):len_trim(parval)).ne."]") &
              call readparams_abort('Arrays must be in square brackets'// &
              ' with elements separated by commas')

         str = trim(adjustl(parval(2:len_trim(parval)-1)))
         
         ! Count how many elements there are
         nparam = 1
         do i=1,len_trim(str),1
            if(str(i:i).eq.",")nparam=nparam+1
         end do
         
         if(nparam.gt.narray) &
              call readparams_abort("Array parameter "//trim(name) &
              //" does not fit in the supplied array!")

         read(str,*,iostat=ios)array(ilower:ilower+nparam-1)
         if(ios.ne.0)call readparams_abort( &
              '  -unable to interpret parameter: '//trim(name))

         if(verbose_params.and.err_params.eq.-1) &
              write(logunit,*)'  -parameter ',trim(name),' = ',array
         if(verbose_params.and.err_params.ne.-1) &
              write(err_params,*)'  -parameter ',trim(name),' = ',array
         
         used_param(jline) = .true.
      else
         nparam = ubound(default,1)-lbound(default,1)+1
      end if

      if(present(nelements))nelements = nparam

      if(write_file)then
         if(nparam.gt.1)then
            write(fmtstring, &
                 '("(a30,"" = "",a1,""["",",i13,"(l13,"",""),l13,""]"")")') &
                 nparam-1
         else
            fmtstring = '(a30,a3,"[ ",l13," ]")'
         end if
         write(out_lun,fmtstring)trim(name)," = ",array(1:nparam)
      end if

      if(make_restartfile)then
         write(restart_lun)nparam
         write(restart_lun)array(1:nparam)
      end if
      if(read_restartfile) then
         read(restart_lun)ntemp
         if(ntemp.ne.nparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         allocate(atemp(1:nparam))
         read(restart_lun)atemp(1:nparam)
         if(any(atemp(1:ntemp).neqv.array(1:nparam)))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         deallocate(atemp)
      end if

      return
    end subroutine get_logical_array_parameter




    subroutine get_string_array_parameter(name,array,default,nelements, &
         required, is_set)

      use log
      implicit none
      character(LEN=*), intent(IN) :: name
      character(LEN=*), intent(OUT), dimension(:) :: array
      integer, optional, intent(OUT) :: nelements
      integer :: iline, jline
      character(LEN=ncharmax) :: str, parname, parval
      integer :: ifirst, ilast, ios
      character(LEN=*), optional, dimension(:) :: default
      integer :: ilower, iupper, narray, nparam
      integer :: i,maxlen,ibegin,iend
      character(LEN=256) :: fmtstring
      logical, optional :: required, is_set
      integer :: ntemp
      character(len=ncharmax), allocatable, dimension(:) :: atemp
    
      if(present(is_set)) is_set = .false.

      if(nlines.lt.1)call readparams_abort('  -parameter file has not been read!')

      ilower = lbound(array,1)
      iupper = ubound(array,1)
      narray = iupper-ilower+1

      jline = -1
      do iline=1,nlines,1
         str = trim(adjustl(line(iline)))
         parname = downcase(trim(str(1:index(str,"=")-1)))
         if(parname.eq.trim(adjustl(downcase(name))))then
            jline = iline
            exit
         end if
      end do

      if(present(required))then
         if(.not.required)then
            if(jline.gt.0)used_param(jline)=.true.
            return
         end if
      end if

      if(jline.lt.0)then
         if(.not.present(DEFAULT)) then
            call readparams_abort( &
                 '  -unable to find parameter '//trim(name)//' in parameter file')
         else
            array = default
            ndefault = ndefault + 1
            write(default_param(ndefault),*)trim(name)//" = ",default 
         end if
      end if

      if(jline.gt.0)then

         if(present(is_set)) is_set = .true.

         str = trim(adjustl(line(jline)))
         ifirst = index(str,"=")+1
         ilast  = index(str,"#")-1 
         if(ilast.eq.-1)ilast=len(str)

         ! Get the parameter value minus any extra spaces at either end
         parval = trim(adjustl(str(ifirst:ilast)))
      
         ! Check that it is enclosed in square brackets
         if(parval(1:1).ne."[".or. &
              parval(len_trim(parval):len_trim(parval)).ne."]") &
              call readparams_abort('Arrays must be in square brackets'// &
              ' with elements separated by commas')

         str = trim(adjustl(parval(2:len_trim(parval)-1)))

         ! Count how many elements there are
         nparam = 1
         do i=1,len_trim(str),1
            if(str(i:i).eq.",")nparam=nparam+1
         end do

         if(nparam.gt.narray) &
              call readparams_abort("  -array parameter "//trim(name) &
              //" does not fit in the supplied array!")
         
         ! Extract elements
         ibegin = 1
         do i = 1, nparam, 1
            iend = ibegin + 1
            do while(str(iend:iend).ne.",".and.iend.lt.len(str))
               iend = iend + 1
            end do
            if(str(iend:iend).eq.",")iend = iend - 1

            array(ilower + i - 1) = trim(adjustl(str(ibegin:iend)))

            ibegin = iend + 2
         end do

         !READ(str,*,iostat=ios)array(ilower:ilower+nparam-1)
         !IF(ios.NE.0)CALL readparams_abort( &
         !     '  -unable to interpret parameter: '//TRIM(name))

         if(verbose_params.and.err_params.eq.-1)then
            write(logunit,*)'  -parameter ',trim(name),' (array of ',nparam,' strings):'
            do i=1,nparam,1
               write(logunit,*)"   "//trim(array(ilower+i-1))
            end do
         end if
         if(verbose_params.and.err_params.ne.-1)then
            write(err_params,*)'  -parameter ',trim(name),' (array of ',nparam,' strings):'
            do i=1,nparam,1
               write(err_params,*)"   "//trim(array(ilower+i-1))
            end do
         end if

         used_param(jline) = .true.
      else
         nparam = ubound(default,1)-lbound(default,1)+1
      end if

      if(present(nelements))nelements = nparam

      if(write_file)then
         write(out_lun,"(a30,a)",advance="no")trim(name)," =  [ "
         do i=1,nparam,1
            write(out_lun,"(a)",advance="no")trim(array(i))
            if(i.lt.nparam)write(out_lun,"(a)",advance="no")" , "
         end do
         write(out_lun,"(a)",advance="yes")" ]"
      end if

      allocate(atemp(1:nparam))
      atemp(1:nparam) = array(1:nparam)
      if(make_restartfile)then
         write(restart_lun)nparam
         write(restart_lun)atemp(1:nparam)
      end if
      if(read_restartfile) then
         read(restart_lun)ntemp
         if(ntemp.ne.nparam)call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
         read(restart_lun)atemp(1:nparam)
         if(any(atemp(1:ntemp).ne.array(1:nparam)))call readparams_abort( &
              '  -restarting but parameter '//trim(name)//' has changed!')
      end if
      deallocate(atemp)

      return
    end subroutine get_string_array_parameter

    
    function downcase(str)
!
! Convert any capitals in an ascii string to lower case.
!
      use log
      implicit none
      character(len=*) :: str
      character :: c
      integer :: i, ic
      character(len=ncharmax) :: downcase
      character(len=ncharmax) :: outstr

      outstr = str

      do i = 1, len_trim(outstr), 1
         c = outstr(i:i)
         if(iachar(c).ge.iachar("A").and.iachar(c).le.iachar("Z"))then
            ic = iachar(c) - iachar("A") + iachar("a")
            outstr(i:i) = achar(ic)
         end if
      end do

      downcase = outstr
      return
    end function downcase


  end module Read_Parameters
